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ABSTRACT 

Dynamic  mechanical  activity  in  a  tunnel  can  be 
measured  as  ground  vibrations  at  offset  distances.  These 
signals  can  be  processed  in  sensing  algorithms  for 
detection,  location,  and  discrimination  of  the  activity.  The 
objective  of  this  work  is  to  demonstrate  that  seismic 
simulations  can  reveal  the  effect  of  the  environment  on 
seismic  energy  as  it  propagates  from  tunnels.  Using 
massively  parallel  high-performance  computers,  the  work 
applies  a  finite-difference  solution  to  the  equations  of 
motion  and  isotropic  stress-velocity  for  viscoelastic 
seismic  propagation.  Results  from  simulations  in  open, 
urban,  and  mountainous  terrain  reveal  the  nature  of 
seismic  waves  as  they  propagate  from  tunnel-digging 
pulses  and  harmonic  sources.  Measures  of  relative  energy 
and  signal  cross  correlation  provide  maps  that  reveal 
locations  of  optimal  sensing.  We  demonstrate  applications 
of  beam  forming  to  monitor  tunnel  activity,  and  conclude 
that  the  simulation  method  produces  realistic  wave-field 
data  for  clarifying  complicated  propagation  phenomena 
and  for  virtual  trials  of  sensing  algorithms. 

1.  INTRODUCTION 

“ The  105th  Military  Police  Battalion  knew  something 
was  amiss.  What  they  discovered  was  breathtaking:  a 
fully  completed  tunnel  that  stretched  357 feet,  longer  than 
a  football  field ’  (Fainaru  and  Shadid,  2005). 

Detection  of  clandestine  tunnels  and  underground 
facilities  is  a  continuing  interest  of  the  US  DoD  (DARPA, 
2006),  Army  (Sabatier  and  Muir,  2006),  and  Customs  and 
Border  Protection  (Rowe,  2006).  Research  continues  in  an 
attempt  to  find  functional  and  reliable  sensing  methods, 
including  seismic  methods. 

Dynamic  activity  in  tunnels  emits  mechanical  energy 
that  propagates  away  in  seismic  waves.  Resulting  ground 
vibrations  can  be  measured  at  offset  distances  and  these 
signals  can  be  used  in  sensing  algorithms  for  detection, 
location,  and  discrimination  of  the  activity.  This  sensing, 
referred  to  as  passive  seismic  detection,  is  complicated  by 
the  topographical,  geological,  infrastructure,  and  noise 
environment  of  the  surrounding  area.  Our  objective  in  this 
paper  is  to  demonstrate  that  high-performance-computing 


(HPC)  simulations  can  reveal  the  effect  of  the 
environment  on  seismic  signatures  from  tunnels,  and  can 
produce  realistic  wave-field  data  for  virtual  trials  of 
sensing  algorithms  and  sensor-placement  decision  tools. 

We  focus  on  results  from  finite-difference  time- 
domain  simulations  of  seismic  wave  propagation  from 
realistic  signature-producing  dynamic  activity  in  tunnels, 
which  we  model  in  open,  urban,  and  mountainous  terrain. 
Using  broadband  sources,  we  characterize  the  frequency 
response  of  our  models,  and  develop  spatial  measures  of 
relative  energy  and  cross  correlation  to  perform  a  virtual 
trial  of  optimal  sensor  placement. 

2.  SEISMIC  ANALYSIS  METHODOLOGY 

The  seismic  simulation  method  is  a  finite-difference 
time-domain  (FDTD)  implementation  of  the  equations  of 
motion  and  isotropic  stress-velocity  equations  for  linear- 
viscoelastic  seismic  propagation.  The  formulation  applies 
stretching  coordinate  transformations,  strain-rate  velocity 
relationships,  first-order  equations  of  motion,  and  linear 
stress-strain  via  stress-rate-strain-rate  equations,  to  derive 
variable-grid  finite-difference  equations  (Ketcham,  et  al., 
2005).  FORTRAN  code  performs  the  computations  on 
HPC  machines  using  four  key  parallelization  features  to 
operate  with  efficient  performance:  equal-domain 
decomposition,  MPI,  efficient  data  exchanges  at  domain 
overlaps,  and  MPIIO. 

We  use  second-order  spatial  finite  differences,  which 
allow  us  to  model  highly  discontinuous  material 
interfaces.  We  can  thus  model  stress-release  at  the  ground 
surface  and  in  structural  members  by  directly  modeling 
air  in  contact  with  these  structures.  This  method  provides 
flexibility  for  complex  topography  and  basic  urban 
settings  within  the  limitations  of  a  rectangular  grid.  The 
variable  grid  in  turn  provides  a  method  to  efficiently 
reduce  a  rectangular  grid’s  staircase-interface  errors.  An 
absorbing  boundary  algorithm  (Cerjan  et  al.,  1985),  used 
in  tandem  with  a  stretched  grid  adjacent  to  the  model 
edges,  minimizes  unwanted  reflections  from  the  edges. 
Seismic  sources  are  introduced  through  body  forces  and 
time  series  detailing  the  force  history. 
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3.  RESULTS 

In  this  section  we  present  results  of  seismic  wave 
propagation  from  repetitive  digging  or  harmonic 
mechanical  sources  inside  tunnels  within  three  models. 
These  models  are  of  flat-layered  open  terrain,  urban 
terrain,  and  mountainous  terrain,  respectively.  For  the 
latter  two  models  we  calculated  the  frequency  response 
function  of  the  ground  surface  to  broadband  sources  in  the 
tunnels,  and  used  the  response  to  examine  locations  for 
sensor  placement.  We  performed  all  of  the  simulations  on 
Cray  XI  computers,  using  multi- streaming  processors  for 
efficient  stress  and  particle-velocity  computations  within 
the  finite-difference  time-stepping  loop. 

3.1.  Propagation  from  tunnel-digging  pulses  in  flat¬ 
layered  terrain 

The  first  model  comprises  a  flat  topographical 
surface,  horizontal  layers  of  soil,  a  short  horizontal  tunnel 
section,  and  a  vertical  tunnel  shaft  open  at  the  surface. 
Both  tunnel  sections  contain  air  in  contact  with  the 
surrounding  ground.  The  soil  layers,  from  the  surface 
downward,  include  1.6  m  of  increasingly  stiff  sandy  soil, 
4.6  m  of  clay,  and  harder  sandy  soil.  We  assigned  the 
viscoelastic  material  properties  and  tunnel  geometry  to 
approximate  near  surface  soil  layers  and  a  constructed 
tunnel  at  an  experimental  site,  where  the  second  author 
has  performed  digging  tests  with  seismic  measurements. 
Fig.  la  illustrates  the  model  tunnel  geometry. 

The  source  in  the  simulation  approximated  force 
pulses  caused  by  a  digger  with  a  hand  instrument  kneeling 
at  the  working  face  of  the  tunnel  (McKenna,  2006; 
McKenna  and  McKenna,  2006).  Vertical  and  horizontal 
force  vectors  combined  to  produce  a  force  into  the  face 
with  a  direction  45  degrees  downward  from  horizontal. 
The  source  time  series  was  a  triangular  pulse  train  with 
0.05-s-duration  pulses  repeating  at  ~1  Hz,  with 
randomness  that  varied  slightly  the  pulse  starting  times 
and  resultant- force  amplitudes.  We  discretized  the  model 
to  produce  accurate  surface  wave  propagation  for  sources 
having  frequency  content  less  than  80  Hz,  and  filtered  the 
input  excitations  so  they  would  not  exceed  this  limit. 

Fig.  lb  contains  a  snapshot  at  t=0.5  s  of  the  simulated 
vertical-particle-velocity  field  on  a  vertical  cross  section 
through  the  tunnel.  The  snapshot  time  is  just  after  the  first 
pulse  application,  and  the  image  reveals  strong  surface 
waves  moving  away  from  the  source  location  and  their 
interference  with  the  shaft.  Signals  extracted  from  the 
simulation  data  reveal  impulsive  vertical  ground-surface 
vibrations  caused  by  the  surface  waves. 

Fig.  2b  plots  simulated  vibrations  from  five  digging 
pulses  superimposed  on  an  ambient  noise  background. 
The  signal  is  from  the  surface  at  8.6  m  from  the  center  of 
the  shaft  in  the  positive  easting  direction.  This  location 
corresponds  to  the  location  of  the  noisy  field 
measurement  shown  in  Fig.  2a,  where  the  digging-pulse 
frequency  was  also  close  to  1  Hz.  The  signals  are  not 


directly  comparable  in  amplitude,  and  their  sampling  rates 
are  different,  but  they  show  qualitatively  comparable 
ground  vibrations  due  to  digging  pulses  that  rise  above 
ambient  background  noise. 


easting  (m) 

Fig.  1.  Flat- terrain  tunnel  model,  (a)  Isometric  view  of 
tunnel  (6.8-m-long  with  lxl-m  cross  section)  and  vertical 
shaft  (6.6-m-deep  with  1.2xl.2-m  cross  section),  (b) 
Zoomed  view  snapshot  of  vertical  particle  velocity,  on 
cross-section  through  tunnel,  caused  by  seismic  waves 
from  digging  pulse.  The  topographic  surface  is  at 
elevation  =  171.75  m.  The  color  bar  is  the  particle- 
velocity  scale.  It  has  a  ±5xl0"8  m/s  range. 

Of  interest  in  real-world  situations  is  the  ability  to 
discern  digging  signatures  not  only  in  an  ambient  noise 
environment,  but  also  in  the  presence  of  specific  noise 
sources.  We  illustrate  the  ability  to  simulate  such 
situations  in  Fig.  3,  which  shows  spectral  content  of  the 
Fig. -2  digging  pulses  when  ground  vibrations  are 
dominated  by  a  nearby  harmonic  source. 

3.2.  Propagation  from  tunnel-digging  pulses  in  urban 
terrain 

The  second  model  is  a  synthesized  urban-terrain 
model,  which  is  described  fully  in  Ketcham  et  al.  (2005). 
This  model  has  two  central  blocks  separated  by  a 
boulevard  with  a  median  and  sidewalks.  Each  central 
block  has  four  five-story  buildings  with  parking  lots  in 
between.  Roads  and  sidewalks  surround  the  central 
blocks.  All  buildings  have  basements.  To  this  model  we 
have  added  an  air-filled  unlined  tunnel  to  represent  a 
clandestine  tunnel  under  construction.  The  tunnel  nearly 
connects  two  buildings. 
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Fig.  2.  Comparison  of  measured  and  simulated  surface 
ground  vibration  from  hand  digging  on  clay  section  of 
tunnel  face  with  -  1-Hz  repetitive  pulses:  (a)  experimental 
and  (b)  synthetic  signals  of  vertical  particle  velocity 
within  noise  background. 
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Fig.  3.  Vertical  ground  vibrations  from  simulated  digging 
pulses  in  the  presence  of  larger  vibrations  from  25-Hz 
generator-induced  noise,  also  from  simulation,  (a) 
Spectrogram  and  (b)  signal  of  the  vibrations  on  surface  35 
m  south  of  the  digging  location.  (The  25-Hz  source  was 
20  m  east  of  the  digging  location).  The  pulses  are  faint  in 
the  spectrogram,  but  evident.  In  the  time  domain  signal 
they  are  not  apparent. 


A  vertical  shaft  from  the  basement  of  one  building 
connects  to  the  horizontal  tunnel,  which  is  9  m  deep  and 
has  a  2x2-m  cross  section.  Figs.  4a  and  4b  illustrate  the 
urban  model  and  soil  layers,  and  Fig.  4c  illustrates  the 
tunnel  from  a  viewpoint  below  the  ground.  The  soil  layers 
have  increasing  stiffness  with  depth,  and  approximate 
subsurface  layering  and  material  properties  at  a  location 
within  Yuma  Proving  Ground.  We  designed  the  force 
vectors  and  time  series  to  approximate  seismic  excitation 
caused  by  digging  on  the  working  face  with  an  air 
hammer,  and  designed  the  model  grid  and  source  filters  to 
provide  accurate  wave  propagation  up  to  48  Hz. 
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Fig.  4.  Urban-terrain  tunnel  model,  (a)  View  from  above, 
(b)  illustration  of  subsurface-soil  layering,  and  (c)  view  of 
tunnel  between  two  buildings,  from  below  the  surface. 
The  two  buildings  with  gray  (rather  than  black) 
highlighting  in  (a)  are  those  shown  in  (c). 

The  excitation  history  was  a  5-s  triangular  pulse-train 
with  0.06-s  pulses  repeating  at  a  mean  frequency  of  33.3 
Hz.  Three  orthogonal  force  vectors  at  the  end  of  the 
tunnel,  with  equal  magnitude  in  the  positive  easting, 
northing,  and  elevation  directions,  combined  to  produce 
the  resultant  forces.  Randomness  in  the  application  time 
and  the  force-resultant  amplitude  approximated  realistic 
variations  of  air-hammer  operation.  The  excitations 
generated  surface  waves  in  the  ground  with  predominant 
wavelengths  slightly  longer  than  10  m.  Fig.  5a  illustrates 
the  vertical  particle  velocity  wave  field  over  the  ground 
surface.  This  snapshot  shows  the  complex  interaction  of 
the  surface  waves  with  the  basements  of  the  buildings, 
and  interactions  with  the  pavement  structures  as  well.  In 
the  soil  just  beneath  the  basements  of  the  buildings,  the 
longer  wavelengths  within  the  deeper  stiffer  soil  are 
evident. 

Figs.  5b  and  5c  plot  the  response  to  the  excitation  at 
-20  m  southwest  from  the  source’s  lateral  position.  Fig. 
5  c  shows  the  time  series  of  vertical  particle  velocity  at 
this  location,  while  Fig.  5b  is  a  spectrogram  that 
illustrates  the  evolving  frequency-domain  power  of  the 
signal.  The  time  series  reflects  the  randomness  given  to 
the  excitation  signal,  while  the  high-amplitude  line  at  -33 
Hz  in  the  spectrogram  shows  that  the  surface  waves  at  the 
signal  location  arrive  at  the  mean  frequency,  and  thus 
with  signature  characteristics,  of  the  source. 
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Fig.  5.  (a)  Zoomed  view  snapshot  of  vertical  particle 
velocity,  on  surface  and  at  base  of  building  foundations, 
from  air-hammer  digging  pulses.  The  particle-velocity 
scale  has  a  ±lxl0"6  m/s  amplitude.  Wave  amplitudes 
saturate  this  scale  near  the  source  but  reveal  the 
propagation  pattern  over  the  domain,  (b)  Spectrogram  and 
(c)  signal  of  vertical  particle  velocity  at  (easting, 
northing)=(197,  120)  m. 

3.3.  Propagation  from  underground-facility  machine 
vibrations  in  mountainous  terrain 

The  third  model  is  a  synthesized  mountainous  setting 
that  includes  an  underground  facility  adjacent  to  a  campus 
of  buildings  that  could  be  apartment  buildings,  industrial 
buildings,  stores,  etc.  Like  the  buildings  in  the  previous 
model,  the  foundations  of  these  buildings  have  sufficient 
depth  (up  to  4  m)  to  affect  surface  wave  propagation.  Fig. 
6a  illustrates  the  model  surface  and  the  locations  of  the 
buildings. 

The  geology  of  the  model  (Ketcham  et  al.,  2004)  has 
the  surface  topography  of  a  mountain  gap.  Its  subsurface 
consists  of  a  soil  basin  at  lower  elevations  and  bedrock 
beneath  a  thin-but-variable-depth  soil  layer  elsewhere. 
The  soil  in  the  basin  extends  irregularly  to  a  depth  as 
much  as  25  m.  Softer  than  the  bedrock,  the  soil  produces 
a  large  contrast  in  amplitude  and  wavelength. 

Fig.  6b  is  a  close-up  view  of  the  campus  and  the 
underground  facility.  The  underground  facility  is  east  of 


the  campus  beneath  an  outcrop  of  bedrock.  It  comprises  a 
central  multi-story  structure  and  a  long  tunnel  with  portals 
emerging  at  the  ground  surface.  The  facility  has  walls  and 
floors  with  material  properties  similar  to  concrete  and 
interior  spaces  with  air  properties.  The  east-west  tunnel  is 
4.5x6x300  m  in  dimension  at  a  600-m  elevation.  The 
central  facility  is  40x40  m  in  plan.  It  has  three  floors, 
each  with  4.5-m  ceilings.  The  lowest  floor  is  at  elevation 
595  m.  Within  this  floor  we  introduced  a  time-varying 
vertical  force  to  approximate  the  excitation  of  a 
harmonically  vibrating  machine  on  a  slab. 


Fig.  6.  Mountainous-terrain  underground-facility  model. 

(a)  View  of  topographic  surface  with  building  locations, 

(b)  close-up  view  of  buildings  and  adjacent  underground 
facility  beneath  outcrop,  and  (c)  zoomed  view  snapshot  of 
vertical  particle  velocity  (m/s)  on  surface  from  machine- 
induced  harmonics.  The  limits  of  the  particle-velocity 
scale  are  ±5xl0"7  m/s. 

We  discretized  the  model  to  analyze  surface  wave 
propagation  for  sources  having  frequency  content  less 
than  25  Hz.  Accordingly,  we  created  the  source  with 
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superimposed  7.5-  and  15-Hz  forces.  These  forces 
resulted  in  a  harmonic  vertical  excitation  of  ±15  kN  at  a 
depth  26  m  below  the  topographic  surface. 

Fig.  6c  illustrates  surface  waves  spreading  from  the 
source  over  the  topography.  The  predominant  features  of 
the  snapshot  are  the  bull’s  eye  pattern  of  propagation,  the 
shortening  of  the  wavelengths  in  the  soil  basin  (to  the  left 
of  the  underground  facility),  and  the  turbulence-inducing 
effect  that  the  building  basements  have  on  the  passing 
waves. 


where  /  =  [  fa ,  fb  ]  designates  the  discrete  frequencies 
from  fa  to  fb  and  *  denotes  the  complex  conjugate.  The 
average  in-band  energy  is: 

fb 

Gf  ( x ,  y)  =  -  V’  /)■ A/  ]  (2) 

Itt  f  ~  fa 

where  A /  is  the  discrete-frequency  difference.  The 
relative  energy,  designated  Gj(x,y) ,  is: 


3.4.  Spatial  processing  of  relative  energy  and 
normalized  cross  correlation 

We  are  currently  researching  the  post-processing  of 
simulated  wave-field  signals  to  calculate  spatial  measures 
of  relative  energy  and  normalized  cross  correlation.  Our 
objective  is  to  map  these  measures  to  view  regions  of 
relatively  high  and  low  signal  amplitude  and  similarity, 
and  to  examine  the  effect  of  complex  propagation 
environments  on  optimal  sensor  placement  for  detection, 
signature,  and  array-processing  analyses.  Here  we  present 
calculations  of  these  measures  from  responses  of  the 
mountainous-  and  urban-terrain  models  to  broadband- 
pulse  excitations.  We  derive  frequency-response  functions 
for  the  surface  nodes,  and  characterize  the  viscoelastic 
systems  excited  previously  by  air-hammer-digging  and 
machine-induced-harmonic  forces,  respectively.  In  effect, 
we  perform  a  trial  of  the  spatial  measures  as  decision- 
support  tools  for  sensor  placement 

Figs.  7a  and  8a  illustrate  a  measure  of  normalized  or 
relative  energy  over  the  surfaces  of  the  mountainous-  and 
urban-terrain  models.  The  measure  is  derived  from  the 
frequency  response  function  between  the  broadband  input 
signal  and  the  output  signal  at  each  node  location  on  the 
model  surface,  which  we  designate  H(x ,  ?/,/),  where  x , 
y  and  /  are  discrete  values  of  easting  and  northing 
coordinates  and  frequency,  respectively.  Our  signal 
processing  code  calculates  the  frequency  response 
function  as  the  ratio  of  the  Fourier  transforms  of  the 
output  signals  to  the  Fourier  transform  of  the  excitation 
signal. 

Using  H(x ,  2/,  /) ,  which  also  equals  the  discrete 
Fourier  transform  of  the  system’s  impulse  response 
function,  the  processing  code  calculates  the  average 
energy  spectral  density  of  the  impulse  response  over 
frequency  bands  of  interest.  For  example,  the  calculation 
in  Fig.  7a  is  for  the  frequency  band  [13.5,  15.7]  Hz,  which 
spans  the  highest-frequency  harmonic  of  the  machine- 
vibration  excitation  in  the  mountainous-terrain  model. 
The  code  calculates  in-band  energy,  average  in-band 
energy,  and  finally  the  relative  energy,  which  we  describe 
here  in  successive  steps.  The  in-band  energy  is: 

G(x,yJ)  =  H*(x,  y,f)H(x,y,  /)  (1) 


Gj(x,  y) 


Gf{x,y) 

median  Gi{x,y)  n 


(3) 


where  the  median  calculation  of  the  denominator  includes 
all  x ,  y  positions  whose  distance  from  the  source 
epicenter  are  within  the  annular  radii  ra  and  rb ;  and  the 
x ,  y  position  of  the  numerator  solution  point  is  within 
this  annulus.  For  Figs.  7a  and  8a  we  chose  A  r  of  the 
radii  to  be  less  than  two  grid-node  spacings.  As  Eq.  3 
indicates,  to  calculate  relative  energy,  the  code  divides  the 
average  in-band  energy  at  each  x,  y  point  by  the  median 
of  all  like  values  within  a  narrow  annulus  from  the  source 
epicenter.  The  annulus  includes  the  point  x,  y  ,  and  so  the 
effect  is  that  distance  from  the  source  is  removed  as  a 
variable. 

Calculating  10  x  log (Gj{x,y))  allows  for  plotting  on 

a  decibel  scale,  with  ±10  dB  spanning  an  order  of 
magnitude.  Using  a  green-yellow-red  color  bar  on  the 
decibel  scale  produces  a  spatial  plot  where,  at  equal 
distance  from  the  source,  green  locations  have  the  highest 
energy  in  the  frequency  band,  i.e.,  the  highest  amplitude, 
and  red  locations  have  the  lowest. 

Figs.  7b  and  8b  illustrate  a  measure  of  the  normalized 
cross  correlation  of  signals  separated  by  a  given  distance, 
respectively,  for  the  mountainous-  and  urban-terrain 
models.  To  derive  the  normalized  cross  correlation  for 
each  node  point,  the  processing  code  first  locates  n  pairs 
of  signals  from  opposite  points  around  a  circle  that  is 
centered  at  the  node  point.  The  2 n  total  signal  locations 
are  equally  spaced  and  thus  separated  by  i t  /  n .  This 
circle  has  diameter  d .  The  processing  provides  a  spatial 
measure  of  similarity  for  signals  separated  by  less  than 
the  distance  d . 

Using  band-limited  signals  within  /  =  [fa,fb],  we 
designate  the  cross  and  auto  correlation  time  series  by 

Rkif(x,y,t),  Rkkf(x,y,t),  and  RUf(x,y,t) ,  where  x  and 

y  define  the  center  of  the  circle  and  k  and  l  designate 
the  two  points  of  a  signal  pair.  The  processing  code 
calculates  the  maximum  of  the  cross  correlation  for  each 
signal  pair  and  divides  this  by  the  product  of  the  auto¬ 
correlation  maxima,  which  occur  at  t  =  0 .  The  average 
of  these  n  calculations  is  the  normalized  cross 

correlation  Rj  (#,  y) ,  i.e., 
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Fig.  7.  Images  of  (a)  relative  energy  (dB)  and  (b) 
normalized  cross  correlation  for  the  mountainous-terrain 
underground-facility  model.  The  calculations  are  for  the 
frequency  band  [13.5,  15.7]  Hz.  The  cross  correlation 
calculations,  at  each  plotted  location,  used  four  signal 
pairs  with  a  spatial  separation  of  120  m. 


Fig.  8.  Images  of  (a)  relative  energy  (dB)  and  (b) 
normalized  cross  correlation  for  the  urban-terrain  tunnel 
model.  The  calculations  are  for  [29.8,  37.3]  Hz.  The 
correlation  calculations  used  four  signal  pairs  at  13.2  m 
separation. 


failing  (m) 

Fig.  9.  Array  configuration  and  locations  (red  circles),  and 
beam  directions  (yellow  lines)  formed  by  beam-forming 
algorithm,  for  array  placements  in  areas  of  (a)  high  and 
(b)  low  relative  energy  and  normalized  cross  correlation, 
plotted  on  a  snapshot  of  vertical  particle  velocity  from  the 
mountainous-terrain  underground-facility  model.  The 
limits  of  the  particle- velocity  scales  are  ±lxl0"7  m/s. 


Fig.  10.  Arrays  (red  circles)  and  beam  directions  (yellow 
lines)  for  arrays  in  areas  of  (a)  high  and  (b)  low  relative 
energy  and  normalized  cross  correlation,  plotted  on  a 
snapshot  of  vertical  particle  velocity  with  noise  from  the 
urban-terrain  tunnel  model.  The  limits  of  the  particle- 
velocity  scales  are  ±lxl0"6  m/s. 
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where  the  points  k  and  l  =  k  +  n  are  sequential 
counterclockwise  points  on  the  circle  each  subtending 
arc-length  i r .  The  values  of  Rj  can  fall  between  0  and  1 . 

We  note  a  related  use  of  normalized  cross  correlation  in 
(Mykkeltveit  et  ah,  1983),  in  which  cross-correlation  data 
were  applied  to  estimate  optimal  array  configurations. 

While  Eq.  4  indicates  time-domain  functions,  our 

calculations  of  Rj  are  in  the  frequency  domain  (Bendat 
and  Piersol,  1986),  beginning  with  H(x,y,f)  extracted 

from  /  =  [fa,  fa].  The  processing  calculates  cross  and 
auto  spectra  of  the  band-limited  impulse  response  at  the 
signal  locations,  and  inverts  to  the  time  domain  to  arrive 
at  the  needed  correlations. 

We  plot  the  result  using  a  green-yellow-red  color-bar 
scale  between  0.65  and  1,  selecting  a  cutoff  below  0.7  to 
indicate  signal  dissimilarity  that  adversely  affects  array 
processing.  Green  locations  have  the  highest  normalized 
cross  correlation  for  the  given  signal  separation  and 
frequency  band,  and  red  locations  have  the  lowest.  Note 
that  processing  with  a  limited  number  of  signal  pairs  at 
each  x ,  y  point  introduces  bias  errors  in  the  plotted 
images,  especially  close  to  the  epicenter.  Fig.  7b  shows 
this  bias  to  have  an  effect  on  normalized-cross-correlation 
values  calculated  with  four  signal  pairs,  producing  a 
radial  pattern  centered  above  the  source  location.  Flat¬ 
layered  models  show  this  effect  to  be  slight,  reducing 

some  values  where  Rj  =  1  to  approximately  0.95. 

Examining  the  plots  of  these  measures  in  Figs.  7  and 
8,  we  can  look  for  regions  with  both  relatively  high 
amplitude  and  normalized  cross  correlation  within  the 
frequency  band  to  support  sensor  locations  for  array 
processing.  For  example,  in  Fig.  7a,  the  effect  of  the 
building  foundations  on  the  amplitude  of  waves 
propagating  from  the  underground  facility  is  evident  in 
the  red-colored  region  to  the  north  and  west  of  the 
building  campus,  which  is  in  the  otherwise  high- 
amplitude  setting  of  the  soil  basin.  Similarly,  this  area  is  a 
region  of  low  normalized  cross  correlation,  as  indicated 
by  Fig.  7b.  Thus  both  the  relative  energy  and  normalized- 
cross-correlation  measures  are  quantitative  indications  of 
the  negative  impact  of  the  building  foundations  on 
seismic  sensing  in  this  region. 

Relative  energy  for  the  urban  tunnel  model  reveals 
the  effect  of  its  buildings  on  propagation  from  the  tunnel. 

The  calculations  spanned  the  frequency  band  /  =  [29.8, 
37.3]  Hz,  which  encompassed  the  air-hammer  pulse 
frequency  of  the  model.  Fig.  8a  shows  the  capacity  of 
building  foundations  to  obstruct  and  direct  energy  in  this 
band,  which  is  evident  in  red-  and  green-colored  regions, 
respectively.  Fig.  8b  shows  the  severe  impact  that 


urbanization,  and  also  the  complexities  of  propagation 
along  the  tunnel,  can  have  on  the  normalized-cross- 
correlation  values.  The  values  are  low  over  most  of  the 
model  surface,  with  the  highest  values  close  to  the  tunnel. 

For  Figs.  7b  and  8b,  we  prepared  the  normalized- 
cross-correlation  maps  using  spatial  separations  of  120 
and  13.2  m,  respectively.  We  chose  these  separations  to 
be  slightly  greater  than  the  wavelengths  of  the  dominant 
waves  in  the  previous  analyses  of  machine  vibrations  and 
tunnel  digging  for  two  reasons:  (1)  our  intended 
application  of  the  maps  was  to  quantify  optimal  array 
positions  for  frequency-wavenumber  (f-k)  beam-forming 
analyses;  and  (2)  the  array  diameter  for  beam-forming 
applications  should  be  at  least  the  length  of  the  longest 
waves  of  interest  (e.g.,  Rost  and  Thomas,  2002).  Using 
Figs.  7  and  8  we  formed  networks  of  beam- forming  arrays 
by  placing  the  arrays  in  “good”  areas,  i.e.,  areas  with 

relatively  high  amplitude  Gj  and  normalized  cross 
correlation  Rj  .  These  locations  are  shown  in  Figs.  9a  and 

10a.  We  also  placed  arrays  in  “bad”  areas  with  relatively 
low  amplitude  and  normalized  cross  correlation.  Figs.  9b 
and  10b  plot  these  arrays.  As  illustrated,  all  arrays 
comprised  eight  locations  on  the  circles  (with  either  120- 
or  13.2-  m  diameter)  plus  the  center  point. 

In  addition  to  the  array  locations,  Figs.  9  and  10  plot 
the  resulting  beams  at  snapshots  in  time  superimposed  on 
respective  images  of  the  wave  field  from  the 
mountainous-terrain  underground-facility  model  and  the 
urban-terrain  tunnel  model.  We  processed  these  beams 
using  a  maximum-slowness  f-k  algorithm  (Hart  and 
Young,  2004).  The  Fig. -9a  arrays,  which  are  from 
locations  with  high  relative  energy  and  normalized  cross 
correlation  in  the  mountainous-terrain  underground- 
facility  model,  show  beams  from  each  array  reliably 
pointing  in  the  direction  of  the  local  wave  arrivals.  The 
beams  also  point  toward  the  source  epicenter,  but  with 
bias  errors  that  derive  from  geology-induced  wave 
refractions.  Without  correcting  for  these  refractions 
(Moran  et  al.,  2001),  source-location  estimates  will  reflect 
these  errors.  Depending  on  the  accuracy  requirement,  a 
source-location  estimate  from  the  Fig. -9a  beams  may  or 
may  not  be  satisfactory.  We  applied  a  median-line¬ 
crossing  algorithm  to  these  and  other  beams  calculated  in 
a  sliding  window  analysis,  which  estimated  the  source 
location  to  within  53-m  accuracy,  or  9  percent  relative 
error  when  normalizing  by  mean  distance  between  the 
source  and  the  arrays.  In  contrast,  the  relative  error  using 
the  Fig. -9b  beams  was  greater  than  75  percent. 

In  our  research,  we  plan  to  use  Gj  and  Rj  with 

measurements  and  estimates  of  noise  power  to  examine 
the  effects  of  complex  geology,  urban  features,  and  noise 
on  sensor  placement  and  performance  as  source  and  noise 
levels  vary  in  time.  As  a  precursor  to  this  work,  we 
superimposed  ambient  noise  onto  the  signals  used  for  the 
beam-forming  analyses  illustrated  in  Figs.  10a  and  10b. 


ma x[Rkif(x,  y,t) \ 

\l Rkk  f  Vi  f  (pi  Vi  0) 


(4) 
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The  signals  are  from  the  simulation  of  air-hammer 
digging  in  the  urban-terrain  tunnel  model.  The  figures 
show  beams  from  arrays  placed  in  “good”  and  “bad” 
sensing  locations,  respectively.  The  wave  field  is  depicted 
with  superimposed  random  ambient  noise.  The  median- 
line-crossing  algorithm,  using  beams  from  the  Fig.- 10a 
arrays,  estimated  the  source  location  to  within  5.8  m 
accuracy  (8  percent  relative  error),  compared  to  9.3  m  (13 
percent  relative  error)  as  the  best  estimate  using  the  Fig.- 
10b  “bad”-location  arrays.  While  this  contrast  in  accuracy 
is  not  as  striking  as  the  previous  example,  both  examples 
show  that  network  source-location  performance  improves 
by  locating  sensors  where  both  amplitude  and  signal- 
similarity  measures  are  relatively  high. 

CONCLUSIONS 

Results  from  seismic  simulations  of  wave 
propagation  from  dynamic  activity  in  tunnels,  which  we 
have  simulated  in  open,  urban,  and  mountainous  terrain, 
reveal  the  pulsating  nature  of  seismic  waves  as  they 
emanate  from  digging  pulses  and  harmonic  vibrations. 
We  find  that  time-domain  characteristics  of  our  simulated 
open-terrain  signals  compare  well  with  field  signals.  In 
mountainous  and  urban  settings,  results  show  the  effect  of 
geology  and  urban  infrastructure  on  propagation  from 
digging-pulse  and  machine-foundation  sources  in  tunnels 
and  underground  facilities.  Our  spatial  post-processing  of 
relative  energy  and  signal  cross  correlation  from  models 
in  mountainous  and  urban  environments  provides  maps  of 
optimal  sensor  performance  in  these  model  domains. 
Application  of  these  maps  is  shown  to  improve  the 
accuracy  of  source  location  by  beam  forming.  We 
conclude  that  the  HPC  seismic  simulation  method 
produces  realistic  wave-field  data,  at  scales  of  interest  for 
DoD  applications,  for  virtual  trials  of  sensing  algorithms 
and  development  of  sensor- array  decision-support  tools. 

ACKNOWLEDGEMENTS 

Funding  support  is  from  DoD  High  Performance 
Computing  Modernization  Program  (HPCMP)  Software 
Applications  Institute  1-01:  Institute  for  Maneuverability 
and  Terrain  Physics  Simulation;  USAERDC  Program 
URBAN;  USA  Rapid  Equipping  Force;  and  Quantum 
Technology  Services,  Inc.,  via  Cooperative  Research  and 
Development  Agreement  “FDTD  Modeling  for  QTSI 
Seismic  Network  Development.”  Computational  support 
is  from  HPCMP  Project  ERDCH1 1304ERD. 

REFERENCES 

Bendat,  J.S.  and  Piersol,  A.G.,  Random  Data,  Analysis 

and  Measurement  Procedures ,  2nd  ed.,  John  Wiley, 

New  York,  1986. 


Cerjan,  C.,  Kosloff,  R.,  and  Reshef,  M.,  A  Nonreflecting 
Boundary  Condition  for  Discrete  Acoustic- Wave  and 
Elastic- Wave  Equations,  Geophys.,  SEG,  Vol.  50, 
1985,  pp.  705-708. 

Defense  Advanced  Research  Project  Agency  (DARPA), 
Passive,  Acoustic,  Seismic  and  Electromagnetic 
Monitoring  (PASEM),  Special  Projects  Office 
Website,  Accessed  June,  2006,  http://www.darpa.mil/ 
spo/programs/pasem.htm. 

Fainaru,  S.  and  Shadid,  A.,  In  Iraq  Jail,  Resistance  Goes 
Underground,  Washington  Post ,  August  24,  2005,  pp. 
A01. 

Hart,  D.  and  Young,  C.,  MatSeis  User's  Manual ,  Ver.  1.9, 
Sandia  National  Laboratories,  Albuquerque,  2004. 

Ketcham,  S.A.,  Moran,  M.L.,  and  Lacombe,  J.,  Seismic 
Waves  from  Light  Trucks  Moving  Over  Terrain, 
Proc.  2004  Users  Group  Conference,  DoD  High 
Performance  Computing  Modernization  Program, 
IEEE  Computer  Society,  pp.  65-70. 

Ketcham,  S.A.,  Lacombe,  J.,  Anderson,  T.S.,  and  Moran, 
M.L.,  Seismic  Propagation  from  Humans  in  Open 
and  Urban  Terrain,  Proc.  2005  Users  Group 
Conference,  DoD  High  Performance  Computing 
Modernization  Program,  IEEE  Computer  Society, 
2005,  pp.  270-277. 

McKenna,  J.R.,  Assessment  of  Seismic- Acoustic 
Techniques  for  Tunnel  Detection,  ERDC  Technical 
Report,  in  review,  76  pp.,  2006. 

McKenna  J.R.  and  McKenna,  M.H.,  Effects  of  Local 
Meteorological  Variability  on  Surface  and 
Subsurface  Seismic-Acoustic  Signals,  in  review, 
Geophysical  Research  Letters ,  2006. 

Moran,  M.L.,  Greenfield,  R.J.,  and  Ketcham,  S.,  Geologic 
Adaptation  for  Seismic  Network  Tracking,  Proc. 
2001  MSS  Battlefield  Acoustic  and  Seismic  Sensing, 
Magnetic  and  Electric  Field  Sensors,  paper  G05, 
Conference  CDROM,  vol.  II,  2002. 

Mykkeltveit,  S.,  Astebol,  K.,  Doornbos,  D.J.,  and 
Husebye,  E.S.,  Seismic  Array  Configuration 
Optimization,  Bull.  Seis.  Soc.  Am .,  SSA,  Vol.  73, 
1983,  pp.  173-186. 

Rost,  S.  and  Thomas,  C.,  Array  Seismology:  Methods  and 
Applications,  Rev.  Geophysics .,  AGU,  Vol.  40,  2002, 
pp.  2-1-2-27. 

Rowe,  R.J.,  Statement  of  Major  General  Richard  J.  Rowe, 
Jr.,  USA,  Director  of  Operations  [J-3,  United  States 
Northern  Command],  Before  the  House  Armed 
Services  Committee,  on  Use  of  the  National  Guard  to 
Support  Border  Security,  US  House  of 
Representatives  Website,  Accessed  June,  2006, 
http://www.house.gov/hasc/schedules/5-24-06Rowe 
Testimony.pdf. 

Sabatier,  J.M.  and  Muir,  T.G.,  Workshop  on  Real-Time 
Detection  of  Clandestine  Shallow  Tunnels,  National 
Center  for  Physical  Acoustics,  Univ.  of  Mississippi, 
NCPA  Report  HB0306-01,  US  Army  Research 
Office,  Grant  No.  W91  INF-06- 1-0001,  April,  2006. 


8 


